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The stochastic approach to the determination of the largest Lyapunov exponent of a many- 
particle system is tested in the so-called mean-field XY-Hamiltonians. In weakly chaotic regimes, 
the stochastic approach relates the Lyapunov exponent to a few statistical properties of the Hessian 
matrix of the interaction, which can be calculated as suitable thermal averages. We have verified 
that there is a satisfactory quantitative agreement between theory and simulations in the disordered 
phases of the XY models, either with attractive or repulsive interactions. Part of the success of the 
theory is due to the possibility of predicting the shape of the required correlation functions, because 
this permits the calculation of correlation times as thermal averages. 
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I. INTRODUCTION 



II. REVIEW OF THE THEORY 



X 



In a recent paper [1] we presented a theory ("the 
stochastic approach" ) that allows to estimate the largest 
Lyapunov exponent of many-particle systems having a 
smooth Hamiltonian. This theory is inspired by and com- 
plementary to the work of Pettini et al. [2,3], and Barnett 
et al. [4]. 

Starting from the linear equations that describe the 
evolution of tangent vectors, we used in [1] the tech- 
niques of stochastic linear differential equations to ob- 
tain an average propagator whose diagonalization pro- 
vides the desired Lyapunov exponent. Even though the 
result obtained in this way is formally (almost) exact, 
approximations must be invoked in practical situations 
to calculate a concrete Lyapunov exponent. 

The aim of this paper is to show that, even if one makes 
several crude approximations, the theory still keeps the 
strength to describe quantitatively some weakly chaotic 
many-particle systems. We do this by comparing the pre- 
dictions of the theory with the outcomings of numerical 
simulations of the Hamiltonian many-particle dynamics. 

The paper is organized in the following way. For rea- 
sons of self-containedness, we begin by presenting a short 
review of the theory in Sect. II. The systems to be stud- 
ied (the mean-field XY-Hamiltonians) are introduced in 
Sect. Ill, where we also work out the predictions of the 
theory for these particular systems. Section IV contains 
the critical comparison of theoretical results with numer- 
ical simulations. We close in Sect. V with some remarks. 



The theory developed in [1] can, in principle, be ap- 
plied to any dynamical system. However, given that the 
approach is perturbative, its success will in general de- 
pend on the choice of the unperturbed system. Here we 
restrict ourselves to perturbations of ballistic motion, i.e., 
Hamiltonians of the type 
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w = 5jX>i+V(9i,...,9Ar), 

i=l 



(1) 



with qi and Pi conjugate position-momentum coordi- 
nates, / the mass (or moment of inertia) of the particles, 
and the interaction V that we assume to be small. 

Let us denote phase space points by x = 
(qi, . . . , qN,pi, ■ ■ ■ ,Pn) and tangent vectors by £ = 
(Sqi, . . . , 5qN, Spi, . . . , Spn) t ■ Differentiating the Hamil- 
ton equations, one obtains the evolution equations for 
tangent vectors: 



(2) 



For a Hamiltonian of the special form (1), the operator 
A has the simple structure 



A(t) 



t/I 
-V(t) 



(3) 



Here V is the Hessian matrix of the potential V, namely 

dqidq-j 
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Once initial conditions xq and £n have been specified, the 
Lyapunov exponent A is found by calculating the limit [5] 



with 



(5) 



We assume that for any initial condition x in phase 
space, the trajectory x(t; x ) is ergodic on its energy shell. 
This implies that A depends only on energy and other sys- 
tem parameters, but not on x , which can then be chosen 
randomly according to the microcanonical distribution. 
There will also be no dependence on initial tangent vec- 
tors, because if £n is also chosen randomly, it will have 
a non-zero component along the most expanding direc- 
tion. Moreover, if the corrections to the exponential law 
implicit in Eq. (5) arc neglected, one can write 



A: 



km -ln(|^;x ,a>)| 2 ) , 

t— »oo Zt 



(6) 



brackets meaning microcanonical averages over xq. This 
approximation will be analyzed in Sect. IV with the help 
of numerical simulations. 

By letting x be a random variable, V(t;x ) becomes 
a stochastic process, and Eq. (2) can be treated as a 
stochastic differential equation. However, as we are in- 
terested in the square of the norm of £, we focus, not on 
Eq. (2) itself, but in the related equation for the evolution 
of the "density matrix" ££ T : 



*Ai(t) = A 1 (t)-(A 1 ) . 



(12) 



A clear exposition of this derivation, together with a de- 
tailed discussion of its domain of validity has been given 
by van Kampen [6] . 

Let i max be the eigenvalue of A which has the largest 
real part. Taking the trace of Eq. (10), one sees that the 
largest Lyapunov exponent A is related to the real part 



of L„ 



by 



A — 2 R- e (^max) 



(13) 



In Eq. (11) we give explicitly only the first two cumu- 
lants, the dots stand for third cumulants and higher order 
ones. If all cumulants were summed up, Eq. (10) would 
be exact in the long-time regime ( » r c [6]. From now 
on, we restrict our analysis to the propagator A trun- 
cated at the second order. The perturbative parameter 
controlling the quality of the truncation is the product 
of two quantities, the "Kubo number" <jt c . The first fac- 
tor, er, characterizes the amplitude of the fluctuations of 
Ai(t). The second, t c , is the correlation time of Ai(t). 

To proceed further one needs the matrix of A in some 
basis. The crudest approximation consists in restricting 
A to the subspacc spanned by the following three matri- 
ces: 



(7) 



the rightmost identity defining the linear superoperator 
A. For the purpose of the perturbative approximations 
to be done, the operator A is split into two parts 



A = Ao + Ai(£) , 



(8) 



where Ao corresponds to the evolution in the absence of 
interactions. In our case A and Ai are associated with 



An = 



1 




and 



Ai 











-V(t) 



> (9) 



respectively (we have set 1=1, but it can be brought 
back at any moment by dimensional considerations). 
Whenever Ai(t) is small, it is possible to manipulate 
Eq. (7) to derive an explicit expression for the evolution 
of the average of ££ T : 



<«'''>(*) = e* A CoCo", 



(10) 



where A is a time-independent superoperator given by 
the perturbative expansion: 

A = Ao + (AO 

+ ^°°dr^A 1 (0e rA ° ( 5A 1 (t-T) e - rAo ) + --- , (11) 



Ii = I J o I , h 




1 



1 

1 



(14) 



This choice is equivalent to a mean-field approximation in 
tangent space [1] . After some lengthy algebra one arrives 
at the corresponding 3x3 matrix for A: 



A = 







-2a 2 r c (3) 



2 

-2fj, 



with the definitions 



a 2 = ^Tr((<5V) 2 ) , 

i k+1) = / drr k f(r), 
Jo 



(15) 

(16) 
(17) 
(18) 



where we have introduced the normalized correlation 
function /(r) 



/(T) = ^T^VW^M) • 



(19) 



[A Lyapunov exponent with the correct units (inverse 
time) is obtained by making the substitutions /i,a — > 
H/I,a/I in Eq. (15).] 
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In the mean-field approximation the Lyapunov expo- 
nent is expressed in terms of the set of four parameters 
H and a 2 Tc k+1 \ k — 0,1,2. The parameters \i and a 
are, respectively, the mean and variance of the stochas- 
tic process V(i), and can in principle be obtained ana- 
lytically by calculating the corresponding microcanonical 
averages. 

The characteristic time ri^ = t c is naturally inter- 
preted as the correlation time of the process V(t). For 
instance, if /(t) is approximately Gaussian, the expan- 
sion of (V(O)V(t)) around r = gives an explicit for- 
mula for the correlation time, namely 



1 




"1 1/2 



TT<7 2 N 



In this case rj 2 ^ and t^' are trivially related to r c : 

T (2) _ ^ 2 
' c _'c ' 



-(3) 



r (3) = l r 3 
7T 



(20) 

(21) 
(22) 



III. THE MEAN-FIELD XY-HAMILTONIAN 

In this section we start the application of the 
perturbative/mean- field theory of Sect. II to a specific 
model. Consider the one-dimensional Hamiltonian 



N N 

H = ^J2 L " + ^H [i -«»(*,-* 



(23) 



This is the so-called mean-field XY-Hamiltonian. It rep- 
resents a lattice of classical spins with infinite-range in- 
teractions. Each spin rotates in a plane and is therefore 
described by an angle < 9i < 2ir, and its conjugate an- 
gular momentum Li, with i — 1,...,N. The constants 
/ and J are the moment of inertia and the interaction 
strength, respectively. (Of course, one can also think of 
point particles of mass I moving on a circle.) 

The mean-field XY-Hamiltonian (HMF) has been ex- 
tensively studied in the last few years (see [7] for a re- 
view). The reasons for the interest in this model are var- 
ious. From a general point of view, the HMF can be con- 
sidered the simplest prototype for complex systems with 
long-range interactions like galaxies and plasmas (in fact, 
the HMF is a descendant of the mass-sheet gravitational 
model [8] ) . But the HMF is also interesting for its anoma- 
lies, be them model-specific or not. Especially worth of 
mention are the long-lived quasi-equilibrium states ob- 
served in the ferromagnetic HMF. These states exhibit 
breakdown of ergodicity, anomalous diffusion, and non- 
Maxwell velocity distributions [9]. The explanation of 
these unusual behaviors may require an extension of the 



standard statistical mechanics, e.g., along the lines pro- 
posed by Tsallis [10]. (Interesting anomalies are also 
present in the anti-ferromagnetic HMF [11]). 

The simplicity of the HMF makes possible a full anal- 
ysis of its statistical properties either in the canonical [8] 
or microcanonical ensemble [12,13]. If interactions are 
attractive (J > 0), the system exhibits a ferromagnetic 
transition at the critical energy E c = 3JN/A. In the case 
J < there is no ordered phase with finite magnetization, 
and, for not too low temperatures, the system behaves 
like the disordered phase of the J > case. However, at 
small energies, some kind of order appears, leading to a 
complex dynamical behavior [11]. 

Lyapunov exponents have been studied with detail, 
both numerically and analytically. For the ferromag- 
netic case, the simulations [7,14,15] show that in the 
magnetized phase the Lyapunov exponent remains finite 
in the thermodynamic limit N — > oo. In contrast, if 
E > E c , A goes to zero when N — > oo; the same behav- 
ior is also observed in the antiferromagnetic case in all 
the energy range [7,16]. The first theoretical studies were 
due to Firpo [17], who derived analytical expressions us- 
ing the geometric method [2,3] (see also the discussion 
in Ref. [16]). Scaling laws in the high-energy regimes 
were derived as well using a random-matrix approach 
[18]. At low temperatures the predictions of the geo- 
metric method for the HMF's Lyapunov exponent are 
not satisfactory. There is solid evidence that the geo- 
metric method predicts wrong scaling laws at low tem- 
peratures for both ferromagnetic and antiferromagnetic 
interactions [14,16]. This means that there are impor- 
tant gaps in the theoretical description of the Lyapunov 
exponent of the HMF (and of many-particle systems, in 
general). Further studies are necessary for understand- 
ing the precise domain of validity of the existing theories: 
this knowledge will be used to make the corresponding 
improvements. 

Going back to the stochastic approach, it was proven 
in [1] that the mean-field approximation is exact in the 
HMF (this agrees with the supersymmetric analysis of 
Tanase- Nicola and Kurchan [19]). In the disordered 
phases of the HMF the fluctuations are small (see be- 
low), so, it is expected that the second-order perturba- 
tive approximation will work well, irrespective of the sign 
of J. 

In order to test the stochastic approach, in its mean- 
field second-order-perturbative version, one has to cal- 
culate the average, variance, and correlation function of 
the Hessian V(i), i.e., Eqs. (16), (17) and (19). These in- 
gredients are microcanonical averages of the appropriate 
observables. In the disordered regimes we will consider, 
microcanonical and canonical averages are equivalent (to 
leading order in N), so, we will prefer the simpler canon- 
ical averaging. Anyway, we will verify numerically that 
canonical and time averages coincide. 
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A. Calculation of /i 

Before embarking in the calculation of canonical aver- 
ages, it is convenient to write the Hamiltonian (23) in 
the simplified form: 



(24) 



i=l 



where we have introduced the magnetization per particle 

JV 



i=i 



with 



fj = (cos 0i, sin 6i) . 



(25) 



(26) 



In terms of m and {r^} the elements of the Hessian ma- 
trix read 



V u = J(m- fi - — ) , 



i +3 



Then one has 



(27) 
(28) 

(29) 



Our next objective is the average (m 2 ). The interacting 
part of the canonical partition function reads 



Z(J3)<x J 



dee 0JNm>/2 _ 



(30) 



The integration over the angles can be reduced to an in- 
tegration over the possible magnetizations: 



Z(0) oc J dmg(m). 



fiJNm 2 j2 



(31) 



Here g(m) represents the density of states in 0-space 
with magnetization m. The problem of finding g(m) has 
a long history and is known as Pearson's random-walk 
problem [20]. There are no simple closed expressions for 
.g(m), but as we are interested in the disordered phases, 
where m is of the order of l/y/W, it suffices to use the 
central-limit approximation: 

g{m) oc c- Nm2 . (32) 

The relative error of this approximation is [20] 

l-^(l-2Nm 2 + 1 -N 2 m ^+... . (33) 

So, we can safely use expression (32) if N » 1 and the 
relevant configurations are such that Nm 2 ~ 1. The last 



condition is synonymous with disorder. It is always satis- 
fied in the antiferromagnetic case. In the case J > it is 
satisfied if the energy is above the critical value e = 3 J/4, 
but not very close to it. 

Putting together the formulas (31) and (32) we get the 
probability distribution of m. To leading order in 1/N 
we have 



P(m) oc exp 

And then we arrive at 

(m 2 )~ 



I ?£) nr 



N (1-/3 J/2) 



(34) 



(35) 



(here and in the following, w means "equal to leading 
order in 1/N"). This expression, together with (24) and 
the equipartition theorem, kT = (Lf /I) , allows us to ob- 
tain the relationship between the temperature T and the 
energy per particle e = E/N: 



kT w 2e - J . 
As a function of e, the average fi reads 

/x/J: 



iV(4e/J-3) ' 



(36) 



(37) 



Combining Eqs. (29) and (35) one sees that \i is a finite- 
temperature correction to the T = oo magnetization. So, 
for fixed N, /i goes to zero as the energy is increased. 
Equation (37) coincides with the result obtained by Firpo 
using a different technique [17]. 



B. Calculation of a 

Given that all degrees of freedom are statistically 
equivalent, the definition of cr 2 , Eq. (17), can be ex- 
pressed as 



a 2 = 



(svy 



ii 



JV 
3=1 



(38) 



the rightmost identity following from the symmetry of V. 
The abovementioned equivalence can be used once more, 
together with the definition of Vij , to obtain 



J2 = ((m-i-i) 2 



1 

+ N 



(m • fi) 

(?i • f 2 ) 2 ) - (?i • r 2 y 



(39) 



To proceed with the evaluation of a 2 one needs the proba- 
bility distributions for two and three particles, ^2(^1, O2) 
and ^3(6*1, 62, O3). Consider first P2, which is just 



P 2 (0i,e 2 ) oc J d6 3 . . . d6 N e P JNm2 / 2 . 



(40) 
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To do the integrations we split the magnetization into 
two parts: the contribution from particles 1 and 2, and 
the remainder: 



m = mi2 + m 3N , 



where 



N 



(?i +r 2 ) , 



m 3A , = 



1 N 



(41) 

(42) 
(43) 



i=3 



We outline the final steps. Like done before, we introduce 
a density of configurations g(m 3 Ar): 



P 2 (0i, 2 ) oc J dm 3N g(m 3N ) exp 



I3JN 



mi2 + m 3 7v) 2 
(44) 



Now invoke the central limit theorem to approximate 
g(m 3 Ar). Switch to polar coordinates. Integrate over 
the polar angle, and then over the modulus m 3 M (the 
upper limit of integration can be extended to infinity). 
The final result is 



P 2 (0 U 6 2 )^A exp[a(f 1 -f 2 )] 

where A is a normalization constant and 

2 

a 



N (4e/J - 3) ' 



(45) 



(46) 



The parameter a is the relevant perturbative quantity in 
this problem. It is small whenever N 3> 1 and the en- 
ergy exceeds the transition value by a finite, large enough, 
amount. (Notice that a = 2/j/J.) 

Proceeding in the same way one also finds the three- 
particle distribution function: 



^3(01, 02, #3) ~ B exp [a(fi • r 2 + r 2 • r 3 + f 3 • fi) 



With the distribution functions P 2 and P 3 in our hands, 
we go back to the calculation of a 2 , Eq. (39). The mo- 
ments of f i • f 2 are immediate: 



(ri • r 2 

\2 



(fl ' f 2 ) 

In addition one has 

(m • fi) ft — - ;r, r_> 



a 

2 ' 
1 

2 ' 



1 

TV 



0(1/JV) . 



Then it is easy to verify that Eq. (39) becomes 
.2 , ! 



J 2 



(m- fi) 



2N 



(48) 
(49) 

(50) 
(51) 



Splitting m into its TV parts, and using that particles are 
statistically equivalent, we arrive at 



11 ~ 1 

J 2 ~ N 



+ ((fi • f 2 ) (fi • f 3 )) . 



(52) 



The three-body average is quickly done (with the help of 
Mathematica [21]) by expanding P 3 in powers of a: 



((fi -f 2 ) (fi -f 3 )) » - . 
So, the final expression for a 2 is 



Na 



^11 + 



(53) 



(54) 



C. Correlation function and correlation time 

Above the transition, the relative importance of the 
interactions decreases with increasing e and N, and the 
dynamics is dominated by the kinetic part of the Hamil- 
tonian. The picture is that of particles rotating almost 
freely during times which are long as compared to the 
mean rotation period. If the system is in equilibrium, 
the dynamics can be modeled by the free-motion equa- 
tions 

e k « e k {t = o) + L fc t/j , i < k < n , (55) 

where {8k(t — 0)} and {Lk} arc independent random 
variables with uniform and Maxwell distributions, re- 
spectively. This is a first approximation valid only during 
short times. Then it is easy to show that the correla- 
tion functions of the elements of the Hessian are directly 
related to the characteristic function of the momentum 
distribution, i.e., 



(47) 2<cos(0 fc -0 J -)Mcos(0 fc -0 j )(O)> « |<exp {-%L h t/I)) | 2 

(56) 



In equilibrium the characteristic function is Gaussian: 



|(exp(-iL fe T/7))| = exp - (r/r*) 



with 



T* 



(57) 



(58) 



After using the definition (18) we arrive at the simplest 
estimate of the correlation time in the high-energy regime 



7T7 



(59) 
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The correlation time t c is of the order of the mean pe- 
riod of rotation. It is independent of the system size 
because it is not directly associated with interactions. 
Of course, interactions are responsible for the Maxwell 
equilibrium distributions, and for the extinction of the 
ballistic regime and its substitution by a random-walk 
one. However, the time scales involved in these processes 
are much longer than r c : the first one is the relaxation 
time for the one-body momentum distribution, the sec- 
ond is the momentum correlation time. Both times grow 
with N, and are not related with r c . 

A more precise estimate for r c can be derived if the 
correlation function is indeed Gaussian. Then the corre- 
lation time is that given by Eq. (20). We will not detail 
the calculation of the canonical averages of Eq. (20) be- 
cause they are very similar to those of Sect. IIIB. Wc 
just show the result: 



h 1 1 


fl + Na/A\ 


/ 4kT 


\1 + Na/8J 



(60) 



very close to the simple estimate of Eq. (59). 

D. The Lyapunov exponent 

Gathering the results of previous sections, one can 
construct the 3x3 matrix A associated to the aver- 
age propagator for the HMF. The Lyapunov exponent 
is extracted from the eigenvalue of A with the largest 
real part. Though the general expression for A might be 
written down explicitly, its content would not justify its 
extension. Notwithstanding, it is worth exhibiting the 
leading term in 1/N. Notice that both /x and a 2 are of 
order 1/N, and that ri fe ' = 0(1) (once the correlation 
function is assumed Gaussian). Then one gets 

A « (^S^j + o (V 2 / 3 ) = O (N- 1 ^) , (61) 

with <7 2 and r c given by Eqs. (54) and (60). The absence 
of \i in the leading-order expression for A is a reflection of 
the fact that in the disordered phases of the HMF, fluc- 
tuations are much larger than the average, i.e. a ^> n, 
and dominate the tangent dynamics. 

IV. NUMERICAL STUDIES 

For comparing the theory with simulations, we will use 
the data existing in the literature [16,22,23] as well as 
some additional data generated by us. We give a succinct 
description of how our simulations were made. Hamil- 
ton's equations of the TV-particle system were evolved 
using the Neri-Yoshida fourth-order symplectic algorithm 
[24]. The time step was fixed through all simulations to 



At = 0.1 (units are such that 1=1 and \J\ = 1). Initial 
conditions were chosen randomly: angles with a uniform 
distribution in [0,27r], and angular momenta Gaussian 
distributed. Then all momenta were shifted by a fixed 
amount to set to zero the total momentum (a constant 
of motion). Finally, velocities were multiplied by an ap- 
propriate factor to fix the total energy at a chosen value 
Ne. The initial one-body distributions generated in this 
way are close to their equilibrium values in the disor- 
dered regimes defined by e > 3J/4 and J > 0, or e > 
and J < 0. Anyway, before doing any "measurements" 
the systems were allowed to relax during a time we call 
t cq , typically t cq e [1000,10000]. Concerning Lyapunov 
exponents, we used Benettin's standard algorithm [5]. 
The initial Euclidian distance between a trajectory and 
its companion were set to c?o = 10~ 6 . We used two al- 
ternative renormalization procedures: (a) The distance 
vector was compressed each time the distance exceeded 
dof; usually with / — 10 or / = 20; (b) Renormaliza- 
tions were made at equally spaced instants, the inter- 
val between successive renormalizations corresponding, 
in average, to expansion factors / w 10, 20. 

Numerically, the Lyapunov exponent is estimated by 
averaging over initial conditions which are propagated 
during a finite time: 

A(t)^l(ln|e(t;x ,£ )| 2 ) . (62) 

We usually considered 10 pairs of randomly-chosen ini- 
tial conditions. Each pair was propagated until a time 
t = tprop, when we judged that (1/t) In |£| had converged 
to a limiting value. Typically, t prop € [1000,5000]. Re- 
member, however, that our theoretical scheme commutes 
logarithm and average. In order to test this approxima- 
tion, we also computed the average 

A*(t)^lln(|^ ;2 ;o,£o)| 2 ) . (63) 

In many cases we run tests to verify that our numer- 
ical results are robust against suitable changes of initial 
conditions, renormalization procedure, or the set of pa- 
rameters {t cq , ip r0 p, At, d ,f}. 



A. Ferromagnetic case 

Figure 1 shows the largest Lyapunov exponent of the 
ferromagnetic HMF as a function of system size N, for 
some selected energies e. The symbols correspond to 
the simulations of Refs. [23] (e = 5.0) and [22] (e = 
1.0,10.0,50.0). We have considered large particle num- 
bers (N > 100) and energies well above the transition 
(e > 1.0) to ensure (i) the validity of the approxima- 
tions invoked in calculating the averages of Sect. Ill, and 
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also (ii) to guarantee that we are in a disordered, quasi- 
ballistic regime. Full lines correspond to the theoret- 
ical Lyapunov exponent obtained by diagonalizing the 
3x3 matrix of Sect. (II), but we could as well have used 
the asymptotic expression of Eq. (61) (minor differences 
would only be visible in the case e = 1). Dotted lines 
correspond to the geometric prediction, as calculated by 
Firpo, i.e., Eqs. (1-3) and (21-22) of Ref. [17]. An inspec- 
tion of Fig. 1 allows us to verify that there is a satisfactory 
agreement between theory and simulation, especially for 
N < 500. For larger systems numerical data deviate 
upward from theory. (Nobre and Tsallis also observed 
deviations from the N~ 1/>3 law in the three-dimensional 
version of the XY-Hamiltonian [25].) This is the opposite 
to natural expectation "the larger the system, the better 
the approximations" . We believe that the reason for this 
deviation may be that the system has not reached micro- 
canonical equilibrium - as far as the Lyapunov exponent 
is concerned. 

It is well known that the HMF relaxes very slowly, e.g., 
Latora, Rapisarda, and Ruffo [14] noticed that for e = 1, 
the system may get trapped in quasi-equilibrium states 
similar to those observed below the transition. A very 
slow relaxation led the same authors to conclude that 
dynamics is ballistic for e = 5.0, N = 1000, t < 10 5 
[26]. If slow relaxation is ruling, the numerical simula- 
tions are likely to provide "quasi-equilibrium" Lyapunov 
exponents, rather than the microcanonical ones. Actu- 
ally, as energies and/or system sizes are increased, the 
system becomes progressively more integrable and devi- 
ations from microcanonical predictions will be necessarily 
observed. 

At this point it is desirable to verify if the canonical cal- 
culations of (j,, er 2 , and r^ k \ reproduce the corresponding 
dynamical averages. In Fig. 2 we show numerical val- 
ues for fi and er 2 , obtained through averaging over 1000 
initial conditions, and in the time window [5000, 10000]. 
One sees that a 2 is very close to the theoretical value, 
within a few percent, but the results for /i show larger 
deviations. We can understand why a 2 is so close to its 
expectation value and fj, is not: The initial value of a 2 
is already very close to the canonical value, within a few 
percent for the cases considered. On the contrary, the 
initial value for fj, is zero, far from its equilibrium value; 
the slow growth of \i is one of the aspects of the slow re- 
laxation to equilibrium. In any case we verified (graphics 
not shown) that the errors in /i and a 2 do not account 
for the deviations of Lyapunov exponents, i.e., feeding 
the theory with the numerical \i and a 2 does not im- 
prove the agreement between theoretical and simulated 
A's. 

In order to discuss the correlation times, let us now 
look at the correlation functions of Fig. 3. Theory and 
simulations agree almost perfectly in the central part of 
the distributions. There are long tails, which cannot be 
appreciated in the figure but might be responsible for 



larger-than-predicted correlation times r^ k \ We recall 
that the quantities are moments of the correlation 
function /(r), and as such, very sensitive to the pre- 
cise shape of the tails We verified that a pure numer- 
ical evaluation of r^ k > fails due to poor convergence of 
the corresponding integrals - this effect is stronger for 
k = 1,2. So, at present we cannot say with certainty if 
the theoretical estimates for r^ k ' , which assume a Gaus- 
sian correlation function, agree with the corresponding 
dynamical ones. To resolve this issue one should make 
a very careful study of the tails of the correlation func- 
tions. Though desirable such analysis exceeds the scope 
of this paper. (Interesting information about correlation 
functions, and, in general, about the geometric method, 
can be found in Refs. [27,28].) 

In passing, let us comment that, for the theory to work, 
/(t) must decay fast enough, e.g., a correlation function 
like r(r) = sinujt/ujt [3,27] would explode the second 
cumulant. 

To conclude the analysis of the ferromagnetic 
possible difference between A and A* [Eqs. (62) and (63)] 
cannot be the explanation for the observed deviations be- 
tween theory and simulations because it is expected that 
A < A* (see below, and Ref. [29]). 

B. Antiferromagnetic case 

Figure 4 shows the results of numerical simulations 
for A and A*. (The numerical results for A, in the 
case N — 100, are consistent with those obtained in 
Ref. [16].) Though not shown, as temperature goes to 
zero (e — > —1/2), the numerical A's go to zero too, the 
empirical scaling law being A cx T 1 / 2 [16] . Notice that the 
relative differences between A and A* are indeed small, 
and decrease with increasing N. We have not analyzed 
the behavior of A* at very low temperatures; in particu- 
lar, we cannot say whether the scaling law A* cx T 1 / 2 is 
verified or not. 

For high energies, i.e., e > 1, our theoretical pre- 
dictions agree reasonably well with simulations. How- 
ever, when the energy tends to its ground-state value, 
s — > —1/2, our theoretical Lyapunov exponent diverges 
dramatically from experiments. These behaviors can be 
easily understood by looking at the energy dependence of 
the perturbative parameter err. Figure (5) shows that, 
for a fixed system size, a 2 remains bounded for all en- 
ergies, growing by a factor of two as T goes from zero 
to infinity. On the other side, r scales with temperature 
like l/Vr, and then, err cx \j\fT. This means that the 
theory must fail as T — > 0. The divergence of the Kubo 
number at T = is a manifestation of the complexity of 
the dynamics close to the antiferromagnetic ground state 
[7,11]. Notwithstanding the overall agreement at high 
temperatures, we must point out a non-resolved source of 
annoyance. At a fixed energy the Kubo parameter scales 
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with the system size like <tt a 1/ >/N. Consequently the 
relative error of the theoretical prediction should decrease 
with increasing N, like 1/y/N. However, e.g., if we take 
e = 7.5, and compare theory versus A*, we see that the 
relative error is about 11%, insensitive to variations of N 
from 100 to 2500. 

For the sake of completeness, we exhibit in Figs. (5) 
and (6) the comparison of canonical and dynamical aver- 
ages for /x, a, and /(r): No significant differences can be 
seen. (Again, the tails of the correlation functions may 
deserve a deeper analysis.) 

Figure (4) also displays the predictions of the geomet- 
ric method [17]. In the high-temperature regime, the rel- 
ative error of the geometric prediction is somewhat larger 
than ours. At low temperatures, the geometric method 
predicts a qualitatively correct behavior, A — > 0, but with 
the wrong scaling law: A oc T 2 [17] instead of A oc T 1 / 2 
[16]. 



V. SUMMARY AND CONCLUDING REMARKS 

In a previous paper [1] we proposed a theoreti- 
cal approach to the determination of the largest Lya- 
punov exponent of many-particle Hamiltonian systems. 
Despite several crude approximations that had to be 
made, it was not completely clear whether the resulting 
perturbative/mean- field script could be carried out for a 
specific system, one of the obstacles being the analytical 
determination of the correlation functions. 

Now we have verified that the stochastic recipe can 
indeed be executed and works satisfactorily in the quasi- 
ballistic regimes of the mean- field XY-Hamiltonians. 
These systems are especially useful for testing the the- 
ory because they make the "mean-field" diagonalization 
exact. However, three additional sources of error have 
to be considered. First of all, we made the simplifica- 
tion of taking the logarithm out of the average [Eq. (6)]. 
That is, using a spin-glass analogy, we are estimating an 
"annealed" Lyapunov exponent A* instead of the usual 
"quenched" A [19]. Quenched averages can be calculated, 
as shown by Tanase-Nicola and Kurchan [19], but they 
require more sophisticated tools, like the replica trick or 
the supersymmetric approach. In the cases considered in 
this paper we verified that the difference between both 
Lyapunov exponents is very small. Then it was not nec- 
essary to go beyond the annealed averaging. 

The second approximation, and most important one, is 
the truncation of the cumulant expansion at the second 
order. By doing so, we introduce a relative error of the or- 
der of the Kubo number ot c . This explains the failure of 
the theory in the low-temperature regime of the antifcr- 
romagnetic XY-Hamiltonian, given that the amplitude of 
the fluctuations a remains bounded but the thermal time 
t c diverges. It is not clear at present if the theory can 



be improved to account for this regime. Further stud- 
ies are necessary. In particular we have to understand 
what kind of dynamical correlations develop as the sys- 
tem approaches the ground state. Another problem that 
requires an answer: We have not observed that the agree- 
ment between theory and simulations betters when the 
size of the system is increased, for fixed energy. This is 
in contradiction with the Kubo number decreasing like 

i/Vn. 

The difficulties of the preceding paragraph may be re- 
lated to the third approximation, which concerns the in- 
teraction representation we have used. In order to isolate 
the fluctuations, one should work in the representation 
associated to the average Hessian. Instead, we have cho- 
sen the free- rotator representation. That is, in construct- 
ing the "free propagator" , we neglected the average Hes- 
sian of the interaction potential. As long as the motion is 
quasi-ballistic, this approximation seems justified. It will 
be certainly wrong, e.g., in the low-temperature phase of 
the ferromagnetic XY, where the quadratic part of the in- 
teraction is really important. We are currently working 
on the implementation of the best interaction represen- 
tation for the infinite-range XY models. The results will 
be presented in a forthcoming paper. 

The mean-field diagonalization is exact for the systems 
considered here. In other cases it represents just a trun- 
cation of the basis for diagonalizing A [Eq. (11)]. Notice, 
however, that we are dealing with a hcrmitian problem, 
i.e., finding the largest eigenvalue of (££ T )(i) [Eq. (10)]. 
Truncation of the basis will produce a lower bound to 
the largest Lyapunov exponent (provided that A is cal- 
culated accurately). So, this problem is analogous to 
finding the ground-state energy of a quantum Hamilto- 
nian. Any small subspace can be considered, the choice 
being guided by the special characteristics of the system 
under study. 

Like its quantum analogue, the Lyapunov problem 
seems reluctant to admit simple general solutions. Each 
class of systems may require special consideration. 
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FIG. 1. Largest Lyapunov exponent of the mean-field 
XY-Hamiltonian (J = 1) as a function of system size 
N, for some selected energies e. From top to bottom, 
e — 1.0,5.0,10.0,50.0, respectively. Symbols correspond to 
numerical simulations of Hamilton's equations. Full lines are 
our theoretical results; dotted lines correspond to the predic- 
tion of the "geometric method" (see text). 
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FIG. 4. Lyapunov exponent of the anti-ferromagnetic 
mean-field Hamiltonian model as a function of energy e. Sym- 
bols are the result of simulations using Benettin's algorithm. 
Open symbols correspond to the usual Lyapunov exponent A 
and hollow symbols to A* (see text). We used t cq — 1000 and 
twin trajectories were propagated during tp rop = 3000. Full 
lines correspond to our theory, dotted ones to the geometric 




method. 
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FIG. 3. Correlation functions for N = 200 and several en- 
ergies, from right to left, e = 1.0,5.0,50.0 (J = 1). Symbols 
correspond to simulations with t cq = 10000 and averaging 
over 100 initial conditions. The lines represent our theoreti- 
cal prediction (Gaussians). 
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FIG. 5. Numerical and analytical values of /i and a 2 for the 
mean-field XY-Hamiltonian model as functions of energy e for 
N = 100 (J = —1). The results shown are time averages in 
the window t e [5000, 10000], and over 100 initial conditions. 
Lines correspond to the corresponding canonical averages. 



FIG. 6. Correlation functions for N = 100 and several en- 
ergies (J = —1). Averaged over 100 initial conditions. Equi- 
libration time, i cq = 10000. Solid lines are our analytical 
results. 
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